// tspcluster - Traveling Salesman Problem
// Copyright (C) 2006  Frederik Carlier
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA

using System;
using System.Collections.Generic;
using System.Text;
using System.Diagnostics;
using System.Collections.ObjectModel;

namespace UGent.Tsp
{
    /// <summary>
    /// Records the transitions made by a Markov Chain. It keeps track of a histogram,
    /// and the transition matrix.
    /// </summary>
    [Serializable()]
    public class Transitions
    {
        private const int MinimumTransitions = 2;

        #region Constructor
        public Transitions()
        {
        }

        public Transitions(int dimension)
            : this(dimension, 0.0, (double)dimension)
        {
        }

        public Transitions(int dimension, double lower, double upper)
            : this(dimension, lower, true, upper, true)
        {
        }

        public Transitions(int dimension, double lower, bool extendLower, double upper, bool extendUpper)
        {
            if (dimension < 0)
                throw new ArgumentOutOfRangeException("dimension", "dimension must be equal or greater than zero.");

            if (dimension > 0 && upper <= lower)
                throw new ArgumentOutOfRangeException("upper", "upper must be greater than lower if dimension is greater than zero.");

            this.dimension = dimension;
            this.normalized = true;

            if (dimension > 0)
            {
                this.lower = lower;
                this.upper = upper;
                this.range = upper - lower;
                this.extendLower = extendLower;
                this.extendUpper = extendUpper;
                this.interval = range / (double)dimension;
                this.transitions = new double[dimension,dimension];
                this.histogram = new long[dimension];
                this.bandWidth = 0;
                this.include = new bool[dimension];
            }
        }
        #endregion

        #region Fields
        private int dimension;
        private double lower;
        private double upper;
        private bool extendLower;
        private bool extendUpper;
        private double range;
        private double interval;
        private long count;
        private int bandWidth; // The number of elements below or above the diagonal that have been set
        private bool[] include;
        private bool normalized; // Normalizing is a one-time only operation.
        private bool merged; // Merging invalidates the histogram.

        private double[,] transitions;
        private long[] histogram;
        #endregion

        #region Properties
        /// <summary>
        /// The number of items recorded.
        /// </summary>
        public long Count
        {
            get
            {
                return count;
            }
        }

        /// <summary>
        /// The dimension (resolution) of the matrix.
        /// </summary>
        public int Dimension
        {
            get { return dimension; }
        }

        /// <summary>
        /// The smallest value in the range that is being recorded.
        /// </summary>
        public double Lower
        {
            get
            {
                return lower;
            }
        }

        /// <summary>
        /// The highest value in the range that is being recorded.
        /// </summary>
        public double Upper
        {
            get
            {
                return upper;
            }
        }

        public bool ExtendLower
        {
            get
            {
                return extendLower;
            }
        }

        public bool ExtendUpper
        {
            get
            {
                return extendUpper;
            }
        }

        /// <summary>
        /// The width of the range that is being recorded.
        /// </summary>
        public double Range
        {
            get
            {
                return range;
            }
        }

        /// <summary>
        /// The spacing between two bins.
        /// </summary>
        public double Interval
        {
            get
            {
                return interval;
            }
        }

        public bool Normalized
        {
            get { return normalized; }
        }

        public bool Merged
        {
            get { return merged; }
        }

        /// <summary>
        /// The histogram.
        /// </summary>
        [CLSCompliant(false)]
        public long[] Histogram
        {
            get
            {
                if (merged)
                    throw new InvalidOperationException();

                return histogram;
            }
        }

        [CLSCompliant(false)]
        public double[,] TransitionMatrix
        {
            get
            {
                return transitions;
            }
        }
        #endregion

        #region General methods
        /// <summary>
        /// Records a new transition
        /// </summary>
        /// <param name="from">The previous value of the Markov Chain.</param>
        /// <param name="to">The newly proposed value.</param>
        /// <param name="type">A value indicating if the move was accepted or not.</param>
        public void AddTransition(double from, double to, TransitionType type)
        {
            if (normalized)
                throw new InvalidOperationException();

            if ((!extendLower && from < lower) || (!extendUpper && from >= upper))
                throw new ArgumentOutOfRangeException("from", "from must be greater than or equal to lower and less than upper");

            if ((!extendLower && to < lower) || (!extendUpper && to >= upper))
                throw new ArgumentOutOfRangeException("to", "to must be greater than or equal to lower and less than upper");

            AddTransition(GetBinIndex(from), GetBinIndex(to), type);
        }

        /// <summary>
        /// Records a new transition
        /// </summary>
        /// <param name="from">The bin index of previous value of the Markov Chain.</param>
        /// <param name="to">The bin index of newly proposed value.</param>
        /// <param name="type">A value indicating if the move was accepted or not.</param>
        public void AddTransition(int from, int to, TransitionType type)
        {
            // Calling this on an empty transition matrix should not fail, but
            // it won't be "normalized" any more - duh.

            if (count > 0 && normalized)
                throw new InvalidOperationException();

            if (dimension == 0)
                return;

            if (normalized)
                normalized = false;

            transitions[from, to] += 1.0;
            // Include[from] and include[to] indicate that a row/column is to be
            // included in the "reduced" matrix (i.e. the matrix to use in the
            // Transition Matrix/Flat Histogram calculations). 
            // It is only set to true if log(transitions[to,from]) > 0 && log(transitions[from,to]) > 0
            // If either is NOT set, but the correspondending reverse transition is positive,
            // we set both
            if (from != to && !(include[from] && include[to]) && (transitions[from, to] >= MinimumTransitions && transitions[from, to] < histogram[from] && transitions[to, from] >= MinimumTransitions && transitions[to, from] < histogram[to]))
            {
                include[from] = true;
                include[to] = true;
            }

            switch (type)
            {
                case TransitionType.Success:
                    histogram[to] += 1;
                    break;
                case TransitionType.Failed:
                    histogram[from] += 1;
                    break;
                default:
                    throw new ArgumentOutOfRangeException("type");
            }

            if (Math.Abs(to - from) > this.bandWidth)
                this.bandWidth = Math.Abs(to - from);

            count++;
        }

        /// <summary>
        /// Gets the bin index for a given value.
        /// </summary>
        /// <param name="value"></param>
        /// <returns></returns>
        public int GetBinIndex(double value)
        {
            if (dimension == 0)
                throw new InvalidOperationException("Cannot get the bin index when the dimension is zero.");

            if ((!extendLower && value < lower) || (!extendUpper && value >= upper))
                throw new ArgumentOutOfRangeException("value", "value must be greater than or equal to lower and less than upper");

            if (value < lower)
                return 0;

            if (value >= upper)
                return dimension - 1;

            // Get the distance between the lower bound of the histogram
            // and the value;
            // In the example: relativeValue = 54 - 10 = 44
            double relativeValue = value - lower;

            // Divide by the interval of the bin, to get the bin number
            // in a real value (ex: 3.14)
            // In the example: binPosition = 44/10 = 4.4
            double binPosition = relativeValue / interval;

            // Use the floor as the real bin number
            // In the example: binNumber = 4
            int binIndex = (int)Math.Floor(binPosition);

            return binIndex;
        }
        #endregion

        #region Histogram methods
        /// <summary>
        /// Returns the percentage of hits that falls into the bin which
        /// contains the specified value.
        /// </summary>
        /// <param name="value">Specifies the bin.</param>
        /// <returns>The percentage of hits in the bin specified by value.</returns>
        public double GetPercentageInBin(double value)
        {
            int binIndex = GetBinIndex(value);

            if (count == 0)
                return 0;

            return histogram[binIndex] / (double)(count);
        }

        /// <summary>
        /// Flattens the histogram.
        /// </summary>
        public void MakeFlatHistogram()
        {
            if (normalized)
                throw new InvalidOperationException();

            for (int i = 0; i < dimension; i++)
                histogram[i] = 1;
        }

        /// <summary>
        /// Clears the histogram.
        /// </summary>
        public void ClearHistogram()
        {
            if (normalized)
                throw new InvalidOperationException();

            for (int i = 0; i < dimension; i++)
                histogram[i] = 0;
        }

        public void Normalize()
        {
            if (normalized)
                return;

            for (int i = 0; i < dimension; i++)
            {
                if (histogram[i] == 0)
                    continue;

                for (int j = 0; j < dimension; j++)
                {
                    transitions[i, j] = transitions[i, j] / histogram[i];
                }
            }

            normalized = true;
        }
        #endregion

        #region Transition matrix methods
        [CLSCompliant(false)]
        public double[] CalculateHistogram()
        {
            return CalculateHistogram(true);
        }

        [CLSCompliant(false)]
        public double[] CalculateHistogram(bool useNative)
        {
            Console.WriteLine("Calculating histogram.");

            if (!normalized)
                throw new InvalidOperationException("not normalized!");

            if (dimension == 0)
                throw new InvalidOperationException("dimension is zero!");

            Console.WriteLine("Banded matrix. Width: {0}.", this.bandWidth);

            // Keeps a list of which rows and columns are included in 
            // the calculation. If a row or column is not included, 
            // the histogram is assumed to be zero.
            Collection<int> reducedIndices = new Collection<int>();
            int reducedDimension = 0;

            for (int i = 0; i < Dimension; i++)
            {
                if (include[i])
                {
                    reducedIndices.Add(i);
                    reducedDimension++;
                }
            }

            if (reducedDimension == 0)
            {
                Console.WriteLine("Reduced dimension is zero!");
                double[] vector = new double[dimension];
                return vector;
            }

            // Create the reduced matrix which we'll ues to calculate.
            // This matrix is a frequency matrix: T(E-E'), normalized to E.
            Console.WriteLine("Creating matrix f");
            double[,] f = new double[reducedDimension, reducedDimension];
            for (int i = 0; i < reducedDimension; i++)
            {
                int nonZeroCount = 0;

                int j_lower = Math.Max(0, i - this.bandWidth);
                int j_upper = Math.Min(i + this.bandWidth + 1, reducedDimension);
                
                for (int j = j_lower; j < j_upper; j++)
                {
                    f[i, j] = transitions[reducedIndices[i], reducedIndices[j]];
                    if (f[i, j] != 0.0)
                        nonZeroCount++;
                }

                if (nonZeroCount < 0)
                    Console.WriteLine(string.Format("WARNING: nonZeroCount for row {0} was {1}", i, nonZeroCount));
            }
            Console.WriteLine("Done!");

            double[] v = new double[reducedDimension];
            double wd, f12, f21, w, sum;
            double[,] cholesky = new double[reducedDimension, reducedDimension];

            for (int bin1 = 0; bin1 < reducedDimension; bin1++)
            {
                v[bin1] = 0;
                wd = 0;

                int bin2_lower = Math.Max(0, bin1 - this.bandWidth);
                int bin2_upper = Math.Min(bin1 + this.bandWidth + 1, reducedDimension);
                
                for (int bin2 = bin2_lower; bin2 < bin2_upper; bin2++)
                {
                    if (bin2 != bin1 && f[bin1, bin2] > 0 && f[bin2, bin1] > 0)
                    {
                        f12 = f[bin1, bin2];
                        f21 = f[bin2, bin1];
                        w = f12 * f21 / (f12 + f21 - 2 * f12 * f21);
                        cholesky[bin2, bin1] = -w;
                        wd = wd + w;
                        v[bin1] = v[bin1] - w * Math.Log(f12 / f21);
                    }
                }

                cholesky[bin1, bin1] = wd;
            }

            for (int bin1 = 0; bin1 < reducedDimension; bin1++)
            {
                wd = cholesky[bin1, bin1];
                if (wd > 0)
                {
                    int upper = Math.Min(bin1 + this.bandWidth + 1, reducedDimension);
                    
                    for (int bin2 = bin1 + 1; bin2 < upper; bin2++)
                    {
                        for (int k = bin1 + 1; k < upper; k++)
                        {
                            cholesky[k, bin2] = cholesky[k, bin2] - cholesky[bin2, bin1] / wd * cholesky[k, bin1];
                        }
                    }

                    for (int k = bin1 + 1; k < upper; k++)
                    {
                        cholesky[k, bin1] = cholesky[k, bin1] / wd;
                    }

                    for (int k = bin1 + 1; k < upper; k++)
                    {
                        cholesky[bin1, k] = cholesky[bin1, k] / wd;
                    }
                }
                else
                {
                    Console.WriteLine("Singulariteit {0} op index {1}", wd, bin1);
                }
            }

            double[] s = new double[reducedDimension];
            // Los het lineair systeem op;
            //    !! solve the linear system
            // DO bin1=1,flathisto_nbins-1
            // s(bin1)=v(bin1) - SUM( s(transition_bin_bounds(1,bin1):bin1-1) &
            // & * cholesky(transition_bin_bounds(1,bin1):bin1-1,bin1) )
            // END DO
            for (int bin1 = 0; bin1 < reducedDimension; bin1++)
            {
                sum = 0.0;
                int lower = Math.Max(0, bin1 - this.bandWidth);
                
                for (int k = lower; k < bin1; k++)
                    sum += s[k] * cholesky[k, bin1];

                s[bin1] = v[bin1] - sum;
            }

            //     DO bin1=1,flathisto_nbins-1
            // IF (cholesky(bin1,bin1)>0) THEN
            // v(bin1)=s(bin1)/cholesky(bin1,bin1)
            // ELSE
            // v(bin1)=0
            // END IF
            // END DO
            for (int bin1 = 0; bin1 < reducedDimension; bin1++)
            {
                if (cholesky[bin1, bin1] > 0)
                    v[bin1] = s[bin1] / cholesky[bin1, bin1];
                else
                    v[bin1] = 0;
            }

            // 
            // s(flathisto_nbins)=0
            // DO bin1=flathisto_nbins-1,1,-1
            // s(bin1)=v(bin1) - SUM( s(bin1+1:transition_bin_bounds(2,bin1)) &
            // & * cholesky(bin1+1:transition_bin_bounds(2,bin1),bin1) )
            // END DO

            s[reducedDimension - 1] = 0;
            for (int bin1 = reducedDimension - 1; bin1 >= 0; bin1--)
            {
                sum = 0.0;
                for (int k = bin1 + 1; k < reducedDimension; k++)
                    sum += s[k] * cholesky[k, bin1];

                s[bin1] = v[bin1] - sum;
            }

            // recenter
            sum = 0.0;
            for(int k = 0; k <reducedDimension; k++)
            {
                sum += s[k];
            }
            for (int k = 0; k < reducedDimension; k++)
                s[k] = s[k] - sum/((double)reducedDimension);

            // Take the exponential
            for (int k = 0; k < reducedDimension; k++)
                s[k] = Math.Exp(s[k]);

            // Normalize
            sum = 0.0;
            for (int k = 0; k < reducedDimension; k++)
                sum += s[k];

            for (int k = 0; k < reducedDimension; k++)
                s[k] = s[k] * (((double)reducedDimension) / sum);

            double[] result = new double[Dimension];
            for (int i = 0; i < reducedDimension; i++)
            {
                result[reducedIndices[i]] = s[i];
            }

            return result;
        }
        #endregion

        public void AddSingleTick()
        {
            double value = 1.0 / ((double)dimension);

            for (int i = 0; i < dimension; i++)
            {
                for (int j = 0; j < dimension; j++)
                {
                    transitions[i, j] = value;
                }
            }
        }

        public Distribution GetDistribution(DistributionEstimator estimator)
        {
            return GetDistribution(estimator, true);
        }
        
        public Distribution GetDistribution(DistributionEstimator estimator, bool normalize)
        {
            Distribution distribution;

            if (estimator == DistributionEstimator.Histogram)
            {
                distribution = new Distribution(lower, extendLower, upper, extendUpper, dimension, (long[])histogram.Clone());
            }
            else
            {
                double[] values = CalculateHistogram(true);
                distribution = new Distribution(lower, extendLower, upper, extendUpper, dimension, (double[])values.Clone());
            }

            if(normalize)
                distribution.Normalize();
                
            return distribution;
        }

        // We can only merge two normalized transition matrices:
        // f^M_{ij} = f^1{ij} + f^2_{ij}
        // Just adding the histograms and the transition matrices
        // doesn't work. We need to work in the *frequency* base.
        // if count = 1, the results of transitions will just be blindly
        // copied over.
        public void Merge(Transitions value)
        {
            if (value == null)
                throw new ArgumentNullException("value");

            // if count = 1, we'll copy the results of transitions into this.
            // Hence, we don't *need* to be normalized!
            if (!Normalized)
                throw new InvalidOperationException("not normalized!");

            if (!value.Normalized)
                throw new ArgumentException("value is not normalized", "transitions");

            if (this.dimension != value.dimension)
                throw new ArgumentOutOfRangeException("value", "Dimensions do not match");

            if (this.lower != value.lower)
                throw new ArgumentOutOfRangeException("value", "Lower bound does not match");

            if (this.upper != value.upper)
                throw new ArgumentOutOfRangeException("value", "upper bound does not match");

            if (this.extendLower != value.extendLower)
                throw new ArgumentOutOfRangeException("value", "extend lower does not match");

            if (this.extendUpper != value.extendUpper)
                throw new ArgumentOutOfRangeException("value", "extendUpper does not match");

            if (this.interval != value.interval)
                throw new ArgumentOutOfRangeException("value", "interval does not match");

            if (this.range != value.range)
                throw new ArgumentOutOfRangeException("value", "range does not match");

            this.count += value.count;

            // If any of the matrices contains a j so that
            // f[i,j]*f[j,i] != null, this matrix (the sum)
            // will satisfy the requirement
            for (int i = 0; i < dimension; i++)
                this.include[i] = this.include[i] || value.include[i];

            this.bandWidth = Math.Max(this.bandWidth, value.bandWidth);
            
            for (int i = 0; i < dimension; i++)
            {
                // If at least one tick has been registered in the new transtion matrix...
                if (value.Histogram[i] > 0)
                {
                    double myCount = this.histogram[i];
                    double otherCount = value.histogram[i];
                    double totalCount = myCount + otherCount;

                    histogram[i] = this.histogram[i] + value.histogram[i];
                    // ... for each row, calculate the average transition rate.
                    
                    int j_lower = Math.Max(0, i-bandWidth);
                    int j_upper = Math.Min(dimension, i+bandWidth+1);
                    
                    for (int j = j_lower; j < j_upper; j++)
                    {
                        transitions[i, j] = myCount / totalCount * this.transitions[i, j] + otherCount / totalCount * value.transitions[i, j];
#if DEBUG
                        Debug.Assert(transitions[i, j] >= 0.0);
                        Debug.Assert(transitions[j, i] <= 1.0);
#endif
                    }
                }
            }

#if DEBUG
            // Do a normalization check
            for (int i = 0; i < dimension; i++)
            {
                double sum = 0.0;
                for (int j = 0; j < dimension; j++)
                    sum += this.transitions[i, j];

                Debug.WriteLine(string.Format("Column sum of column [{0}]: {1} (should be 0 or 1)", i, sum));
            }
#endif

            // However, it is possible that a non-zero f[i,j] in one matrix
            // is matched by a non-zero f[j,i] in another. Hence, we'll need 
            // to re-evalute
            for (int i = 0; i < dimension; i++)
            {
                // If we already know that the requirement is satisfied, 
                // continue;
                if (this.include[i])
                    continue;

                // If there are no events, we can skip this one, too!
                if (this.histogram[i] == 0)
                    continue;

                int j_min = Math.Max(0, i - this.bandWidth);
                int j_max = Math.Min(dimension - 1, i + this.bandWidth);

                for (int j = j_min; j <= j_max; j++)
                {
                    if (i != j && this.transitions[i, j] > 0.0 && this.transitions[j, i] > 0.0 && this.transitions[i, j] < 1.0 && this.transitions[j, i] < 1.0)
                    {
                        this.include[i] = true;
                        this.include[j] = true;
                    }
                }
            }
        }
        
        public Transitions Clone()
        {
            Transitions value = new Transitions(this.dimension, this.lower, this.extendLower, this.upper, this.extendUpper);
            value.bandWidth = this.bandWidth;
            value.count = this.count;
            value.histogram = (long[])this.histogram.Clone();
            value.include = (bool[])this.include.Clone();
            value.interval = this.interval;
            value.merged = this.merged;
            value.normalized = this.normalized;
            value.range = this.range;
            value.transitions = (double[,])this.transitions.Clone();
            
            return value;
        }
    }
}
